function [ dat wnorms wnorm wnorms_nai wnorm_nai ] = get_voxel_recon_timecourse( S )

% [ dat wnorms wnorm wnorms_nai wnorm_nai ] = get_voxel_recon_timecourse( S )
%
% MWW 2012

sensor_data=S.sensor_data;
source_recon_results=S.source_recon_results;
chanind=S.chanind;
class_samples_inds=S.class_samples_inds;
voxind=S.voxind; %=first_level_results.mask_indices_in_source_recon(indind)

if ~strcmp(source_recon_results.recon_method,'none'), % not sensor space analysis
    if isfield(source_recon_results.BF.inverse.W,'MEG') % added by DM
        NK=length(source_recon_results.BF.inverse.W.MEG);
    else
        NK=length(source_recon_results.BF.inverse.W.EEG);
    end
else
    NK=1;
end;

Ntrials=size(class_samples_inds{1},3);
Ntpts=length(find(source_recon_results.samples2use));

dat=nan(Ntrials,Ntpts);
wnorms=nan(Ntrials,Ntpts);
wnorms_nai=nan(Ntrials,Ntpts);
wnorm=zeros(NK,1);
wnorm_nai=zeros(NK,1);
weights=zeros(NK,length(chanind));

for kk=1:NK,  

    if strcmp(source_recon_results.recon_method,'none'), % sensor space analysis

        % for sensor space - just compute for all sensors
        weights(kk,voxind)=1;
    else
        
        if isfield(source_recon_results.BF.inverse.W,'MEG')  % added by DM
            weights(kk,:)=source_recon_results.BF.inverse.W.MEG{kk}{voxind};
        else
            weights(kk,:)=source_recon_results.BF.inverse.W.EEG{kk}{voxind};
        end
    end;
    
    if sum(isnan(squash(weights(kk,:))))==0,

        if ~strcmp(source_recon_results.recon_method,'none'),
            % this one represents the uncertainty and will be
            % applied equally to the data and regressors, is
            % therefore irrelevant if NK=1
            
            if isfield(source_recon_results.BF.inverse.W,'MEG')  % added by DM
                cov=source_recon_results.BF.features.C.MEG{kk};
            else
                cov=source_recon_results.BF.features.C.EEG{kk};
            end
            
            wnorm(kk)=weights(kk,:)*cov*weights(kk,:)';
            % this one represents a scaling of just the data
            wnorm_nai(kk)=weights(kk,:)*weights(kk,:)';   

        else
            wnorm(kk)=1;
            wnorm_nai(kk)=1;
        end;
    else,
        disp('NAN weights');
    end;
    
end;

if(S.use_class_probs),
    class_prs=permute(S.class_prs,[2 3 1]);
    dims=size(class_prs);
    class_prs=reshape(class_prs,[dims(1)*dims(2) dims(3)]);
    weights2=class_prs*weights;
    
    sensor_data2=permute(sensor_data,[2 3 1]);
    dims=size(sensor_data2);
    sensor_data2=reshape(sensor_data2,[dims(1)*dims(2) dims(3)]);
    
    dat=zeros(size(weights2,1),1);
    wnorms=zeros(size(weights2,1),1);
    wnorms_nai=zeros(size(weights2,1),1);
    for ii=1:size(weights2,1),
        dat(ii)=weights2(ii,:)*sensor_data2(ii,:)';
        if isfield(source_recon_results.BF.inverse.W,'MEG')  % added by DM
            cov=source_recon_results.BF.features.C.MEG{kk};
        else
            cov=source_recon_results.BF.features.C.EEG{kk};
        end
        wnorms(ii)=weights2(ii,:)*cov*weights2(ii,:)';
        wnorms_nai(ii)=weights2(ii,:)*weights2(ii,:)';
    end;
    dat=reshape(dat,dims(1),dims(2))';
    wnorms=reshape(wnorms,dims(1),dims(2))';
    wnorms_nai=reshape(wnorms_nai,dims(1),dims(2))';
else,

    for kk=1:NK,
        for tri=1:Ntrials, % indexes trials
            timeinds=find(class_samples_inds{kk}(1,:, tri)); % time indices for class kk
            if ~isempty(timeinds),
                dat(tri,timeinds)=weights(kk,:)*sensor_data(:,timeinds,tri);
                %wnorms(tri,timeinds)=wnorm(kk);    
                wnorms(tri,timeinds)=wnorm(kk);                   
                wnorms_nai(tri,timeinds)=wnorm_nai(kk); 
            end;
        end;
    end;
end;

end

